home *** CD-ROM | disk | FTP | other *** search
/ ADA Programming Guide / ADA Programming Guide.iso / ada_gwu / farith.c < prev    next >
C/C++ Source or Header  |  1996-01-30  |  13KB  |  509 lines
  1. /*
  2.  * Copyright (C) 1985-1992  New York University
  3.  * 
  4.  * This file is part of the Ada/Ed-C system.  See the Ada/Ed README file for
  5.  * warranty (none) and distribution info and also the GNU General Public
  6.  * License for more details.
  7.  
  8.  */
  9. /*
  10.  * multi-precision multiplication, division and exponentiation
  11.  * for fixed point computations.
  12.  */
  13.  
  14. #include <stdlib.h>
  15. #include "config.h"
  16. #include "int.h"
  17. #include "ivars.h"
  18. #include "intbp.h"
  19. #include "miscp.h"
  20. #include "farithp.h"
  21.  
  22. #define int_copy(u,v)   for(i=0;i<=FIX_LENGTH;i++) v[i] = u[i]
  23. #define int_conv(v,w)   w[0]=1, w[1]=v
  24.  
  25. static void int_norm(int *);
  26. static void int_print(int *);
  27.  
  28. #ifdef DEBUG_FARITH
  29. void int_print();
  30. #endif
  31.  
  32. /*
  33.  * M U L T I P L E     P R E C I S I O N     I N T E G E R
  34.  *
  35.  *       A R I T H M E T I C       P A C K A G E
  36.  *
  37.  *             Robert B. K. Dewar
  38.  *              June 16th, 1980
  39.  *
  40.  * This package of routines implements multiple precision integer
  41.  * arithmetic using what are called the "classical algorithms" in
  42.  * Knuth "Art of Programming", Volume 2, Section 4.3.1. The style
  43.  * of the algorithms follows Knuth fairly closely, and section
  44.  * 4.3.1 can be consulted for further theoretical details.
  45.  *
  46.  * Multiple precision integers are represented as tuples of small
  47.  * integers in the range 0 to FIX_BAS - 1, where FIX_BAS is a power of 2
  48.  * (actually 2 ** FIX_DIGS) which is the base used to represent the
  49.  * long integers. Essentially the successive elements of the tuple
  50.  * are the digits of the representation in base FIX_BAS. All integers
  51.  * are normalized so that the first digit is non-zero, except in
  52.  * the case of zero itself. All values are positive.
  53.  *
  54.  
  55.  * The constants FIX_BAS and FIX_DIGS must be defined as global constants
  56.  * in a program using these routines. It is assumed that the value
  57.  * FIX_BAS * FIX_BAS - 1 can be represented as an integer value.
  58.  
  59.  * The following routines are provided:
  60.  
  61.  *      int_norm        integer normalization
  62.  *    int_div        integer division
  63.  *    int_mul        integer multiplication
  64.  *    int_tom        conversion to multi-precision
  65.  *    int_tol        conversion to long integer
  66.  */
  67.  
  68.  
  69. void int_tom(int *v, long n)                                    /*;int_tom*/
  70. {
  71.     /* Convert a positive integer to a multiple precision integer. */
  72.  
  73.     int        d;
  74.  
  75.     /* Build up v in groups of FIX_BAS digits until n is depleted. */
  76.  
  77.     d = v[0] = 5;
  78.     while(d) {
  79.         v[d--] = n % FIX_BAS;
  80.         n /= FIX_BAS;
  81.     }
  82.     int_norm(v);
  83. }
  84.  
  85. void int_mp2(int *u, int p)    /*;int_mp2*/
  86. {
  87.     /* Multiplication by a power of 2. (Shift) */
  88.  
  89.     int     s;
  90.     int     m;
  91.     int     i,f,t,c;
  92.  
  93.     m = p % FIX_DIGS;
  94.     s = p / FIX_DIGS;
  95.     c = 0;
  96.  
  97.     if (m) { /* needs multiplication */
  98.         f = 2;
  99.         for (i=1;i<m;i++) f *= 2; /* f = 2 ** m */
  100.         if (f * u[1] > FIX_BAS) { /* carry: needs extension */
  101.             for (i = u[0]; i > 0; i--) {
  102.                 t = u[i] * f + c;
  103.                 u[i+1] = t % FIX_BAS;
  104.                 c = t / FIX_BAS;
  105.             }
  106.             u[1] = c;
  107.             u[0] += 1;
  108.         }
  109.         else {
  110.             for (i = u[0]; i > 0; i--) {
  111.                 t = u[i] * f + c;
  112.                 u[i] = t % FIX_BAS;
  113.                 c = t / FIX_BAS;
  114.             }
  115.             /* carry is zero now */
  116.         }
  117.     }
  118.     /* now we just have to add some zeros at the end */
  119.     if (s) { /* needs shift */
  120.         for (i=1;i<=s;i++)
  121.             u[u[0]+i] = 0;
  122.         u[0] += s;
  123.     }
  124. }
  125.  
  126. void int_mul(int *u, int *v, int *w)                                /*;int_mul*/
  127. {
  128.     /* Multiply unsigned integers, cf Knuth 4.3.1 Algorithm M
  129.      *
  130.      * Multiplication is similar to, but not identical with, the
  131.      * usual pencil and paper algorithm. The main difference is
  132.      * that the sum is accumulated as we go along, rather than
  133.      * forming all the partial sums and adding them up at the end.
  134.      *
  135.      */
  136.  
  137.     int        i, j, k, t;
  138.  
  139.     /* Initialize result to all zeroes(actually it is only absolutely
  140.      * necessary to initialize the last #v digits of w to zero).
  141.      */
  142.  
  143. #ifdef DEBUG_FARITH
  144.     printf("int_mul\n"); 
  145.     int_print(u); 
  146.     int_print(v);
  147. #endif
  148.     w[0] = u[0] + v[0];
  149.     for (i = 1; i <= FIX_LENGTH; i++) w[i] = 0;
  150.  
  151.     /* Outer loop, through digits of multiplier */
  152.  
  153.     for (j = v[0]; j > 0; j--) {
  154.         /* The inner loop, through the digits of the multiplicand, is
  155.          * similar to the inner loop of the addition, except that the
  156.          * product is added in, and the carry, k, can exceed 1.
  157.          */
  158.  
  159.         k = 0;
  160.         for (i = u[0]; i > 0; i--) {
  161.             t = u[i] * v[j] + w[i + j] + k;
  162.             w[i + j] = t % FIX_BAS;
  163.             k = t / FIX_BAS;
  164.         }
  165.  
  166.         /* The final step in the inner loop is to store the final
  167.          * carry in the next position in the result.
  168.          */
  169.  
  170.         w[j] = k;
  171.     }
  172.  
  173.     /* The result must be normalized(there could be one leading zero),
  174.      * and then the result sign attached to the returned value.
  175.      */
  176.  
  177.     int_norm(w);
  178. #ifdef DEBUG_FARITH
  179.     int_print(w);
  180. #endif
  181. }
  182.  
  183. long int_tol(int *t)                                            /*;int_tol*/
  184. {
  185.     /* Convert a multiple precision integer to a C long integer, if possible.
  186.      * Otherwise, return 'OVERFLOW'.
  187.      *
  188.      * First check if it overflows.
  189.      */
  190.  
  191.     long    x;
  192.     int        i;
  193.  
  194.     arith_overflow = 0;        /* reset overflow flag */
  195.  
  196.     if (t[0] > 5) {
  197.         arith_overflow = 1;
  198.         return 0;
  199.     }
  200.     if (t[0] == 5 && t[1] >= 8) { /* t >= 8 * 2**28 = 2**31 */
  201.         arith_overflow = 1;
  202.         return 0;
  203.     }
  204.  
  205.     x = t[1]; /* set first part of result */
  206.     for (i = 2; i <= t[0]; i++)
  207.         x = FIX_BAS * x + t[i];
  208.  
  209.     return (x);
  210. }
  211.  
  212. void int_div(int *u, int *v, int *q)                            /*;int_div*/
  213. {
  214.     /* Obtain quotient and remainder of signed integers,
  215.      * cf Knuth 4.3.1 Algorithm D
  216.      * Result is returned as a tuple [quotient,remainder].
  217.      *
  218.      * This is by far the most difficult of the four basic operations.
  219.      * This is because the paper and pencil algorithm involves certain
  220.      * amounts of guess work which cannot be programmed directly. The
  221.      * approach(which is analyzed in detail in Knuth) is to reduce
  222.      * the guess work by computing a rather good guess at each digit
  223.      * of the result, and then correcting if the guess is wrong.
  224.      *
  225.      * If the divisor is zero, return om.
  226.      */
  227.  
  228.     int        i, j, k, du, p, c, d;
  229.     int        rr, n, m, qe;
  230.  
  231.     /* A special case, if u is shorter than v, the result is 0 */
  232.  
  233.     if (u[0] < v[0]) {
  234.         int_conv(0,q);
  235.         return;
  236.     }
  237.  
  238.     /* The case of a one digit divisor is treated specially. Not only
  239.      * is this more efficient, but the general algorithm assumes that
  240.      * the divisor contains at least two digits.
  241.      */
  242.  
  243.     if (v[0] == 1) {
  244.         q[0] = u[0];
  245.  
  246.         /* Basically this case is straight-forward. Since we can represent
  247.          * numbers up to FIX_BAS * FIX_BAS - 1, we can do the steps of the
  248.          * division exactly without any need for guess work. The division is
  249.          *  then * done left to right (essentially it is the short division
  250.          * case).
  251.          */
  252.         rr = 0;
  253.         for (j = 1; j <= u[0]; j++) {
  254.             du = rr * FIX_BAS + u[j];
  255.             q[j] = du / v[1];
  256.             rr = du % v[1];
  257.         }
  258.     }
  259.     /* Here is where we must do the full long division algorithm */
  260.  
  261.     else {
  262.         n = v[0];
  263.         m = u[0] - v[0];
  264.         u[0] += 1; /* extend u */
  265.         for(i=u[0];i>1;i--) u[i] = u[i-1];
  266.         u[1] = 0;
  267.         q[0] = m+1;
  268.  
  269.         /* The first step is to multiply both the divisor and dividend
  270.          * by a scale factor. Obviously such scaling does not affect
  271.          * the quotient. The purpose of this scaling is to ensure that
  272.          * the first digit of the divisor is at least FIX_BAS / 2. This
  273.          * condition is required for proper operation of the quotient
  274.          * estimation algorithm used in the division loop. Note that
  275.          * we added an extra digit at the front of the dividend above.
  276.          */
  277.  
  278.         d = FIX_BAS /(v[1] + 1);
  279.  
  280.         c = 0;
  281.         for (j = u[0]; j > 0; j--) {
  282.             p = u[j] * d + c;
  283.             u[j] = p % FIX_BAS;
  284.             c = p / FIX_BAS;
  285.         }
  286.  
  287.         c = 0;
  288.         for (j = v[0]; j > 0; j--) {
  289.             p = v[j] * d + c;
  290.             v[j] = p % FIX_BAS;
  291.             c = p / FIX_BAS;
  292.         }
  293.  
  294.         /* This is the major loop, corresponding to long division steps */
  295.  
  296.         for (j = 1; j <=(m + 1); j++) {
  297.             /* Guess the next quotient digit by doing a division based on the
  298.              * leading digits. This estimate is never low and at most 2 high.
  299.              */
  300.             qe =(u[j] != v[1])
  301.               ? ((u[j] * FIX_BAS + u[j + 1]) / v[1]) :(FIX_BAS - 1);
  302.  
  303.             /* The following loop refines this guess so that it is almost always
  304.              * correct and is at worst one too high(see Knuth for proofs).
  305.              */
  306.  
  307.             while((v[2] * qe) >
  308.               ((u[j] * FIX_BAS + u[j + 1] - qe * v[1]) * FIX_BAS + u[j + 2])) {
  309.                 qe -= 1;
  310.             }
  311.  
  312.             /* Now(for the moment accepting the estimate as correct), we
  313.              * subtract the appropriate multiple of the divisor. This is
  314.              * similar to the inner loop of the multiplication routine.
  315.              */
  316.             c = 0;
  317.             for (k = n; k > 0; k--) {
  318.                 du = u[j + k] - qe * v[k] + c;
  319.                 u[j + k] = du % FIX_BAS;
  320.                 c = du / FIX_BAS;
  321.                 if (u[j + k] < 0) {
  322.                     u[j + k] += FIX_BAS;
  323.                     c -= 1;
  324.                 }
  325.             }
  326.             u[j] += c;
  327.  
  328.             /* If the estimate was one off, then u(j) went negative when the
  329.              * final carry was added above. In this case, we add back the
  330.              * divisor once, and adjust the quotient digit.
  331.              */
  332.  
  333.             if (u[j] < 0) {
  334.                 qe -= 1;
  335.  
  336.                 c = 0;
  337.                 for (k = n; k > 0; k--) {
  338.                     u[j + k] += v[k] + c;
  339.                     if (u[j + k] >= FIX_BAS) {
  340.                         u[j + k] -= FIX_BAS;
  341.                         c = 1;
  342.                     }
  343.                     else
  344.                         c = 0;
  345.                 }
  346.  
  347.                 u[j] += c;
  348.             }
  349.  
  350.             /* Store the next quotient digit */
  351.  
  352.             q[j] = qe;
  353.         }
  354.     }
  355.  
  356.     /* All done, except for normalizing the quotient
  357.      * to remove leading zeroes and supplying the
  358.      * proper result sign to the returned values.
  359.      */
  360.  
  361.     int_norm(q);
  362. }
  363.  
  364. static void int_norm(int *u)        /*;int__norm*/
  365. {
  366.     /* Remove leading zeroes from calculated value
  367.      *
  368.      * The representation used in this package requires that all integer
  369.      * values be normalized, i.e. the first digit of any stored value
  370.      * must be non-zero except for the special case of zero itself. The
  371.      * normalize routine is called to ensure this condition is met.
  372.      *
  373.      */
  374.  
  375.     int        i, j;
  376.  
  377.     if (u[0] == 1 || u[1] != 0)
  378.         return;
  379.  
  380.     for (i = 2; i <= u[0]; i++) {
  381.         if (u[i]) {
  382.             u[0] = u[0] -(i - 1);
  383.             for (j = 1; j <= u[0]; j++)
  384.                 u[j] = u[i++];
  385.             return;
  386.         }
  387.     }
  388.     /*  Return zero if all components zero */
  389.     u[0] = 1;
  390.     return;
  391. }
  392.  
  393. /* debugging and test procedures */
  394. static void int_print(int *u)                                    /*;int_print*/
  395. {
  396.     /* Dump multi-precision integer to standard output */
  397.  
  398.     int        i,n;
  399.  
  400.     n = u[0];
  401.     if (n <= 0)
  402.         chaos("ill-formatted arbitrary precision integer - lengt <=0");
  403.     printf("integer: %d components\n", u[0]);
  404.     printf("%3d %*d\n", 1, DIGS+1, u[1]); /* allow for possible sign */
  405.     for (i = 2; i <= u[0]; i++)
  406.         printf("%3d  %0*d\n", i, DIGS, u[i]);
  407. }
  408.  
  409. #ifdef DEBUG_FARITH
  410. #define NUMBER 20
  411. #define LENGTH 10
  412. main ()
  413. {
  414.     int pow5[100];
  415.     int i,j,c,v;
  416.     pow5[0] = 1;
  417.     pow5[1] = 1;
  418.     printf("static int pow5[%d] [%d] = {\n",NUMBER+1,LENGTH+1);
  419.     for (i=0;i<NUMBER;i++) {
  420.         printf("    {");
  421.         for (j=0;j<=pow5[0];j++) {
  422.             printf(" %3d,",pow5[j]);
  423.         }
  424.         for (;j<LENGTH;j++) {
  425.             printf("   0,");
  426.         }
  427.         printf("   0 },\n");
  428.         for (c=0,j=pow5[0];j;j--) {
  429.             v = 5 * pow5[j] + c;
  430.             pow5[j+1] = v % FIX_BAS;
  431.             c = v / FIX_BAS;
  432.         }
  433.         if (c) {
  434.             pow5[1] = c;
  435.             pow5[0] = pow5[0] + 1;
  436.         }
  437.         else {
  438.             for (j=1;j<=pow5[0];j++) {
  439.                 pow5[j] = pow5[j+1];
  440.             }
  441.         }
  442.     }
  443.     printf("    {");
  444.     for (j=0;j<=pow5[0];j++) {
  445.         printf(" %3d,",pow5[j]);
  446.     }
  447.     for (;j<LENGTH;j++) {
  448.         printf("   0,");
  449.     }
  450.     printf("   0 }\n};\n");
  451.     if (pow5[0] >= LENGTH) {
  452.         printf("Fatal: LENGTH is too short\n");
  453.         exit();
  454.     }
  455.  
  456.     printf("pow_of5(v,p)\t\t\t/*;pow_of5*/\n");
  457.     printf("/* This procedure is generated automatically \n");
  458.     printf(" * and therefore should not be modified.\n */\n");
  459.     printf("int *v,p;\n{\nint i,*u;\n");
  460.     printf("\tif (p < 0 || p > %d) {\n",NUMBER);
  461.     printf("\t\tv[0] = v[1] = 1;\n");
  462.     printf("\t\traise(SYSTEM_ERROR,\"power of 5 too large\");\n\t}\n");
  463.     printf("\tu = pow5[p];\n");
  464.     printf("\tfor (i=0;i<=u[0];i++) v[i] = u[i];\n}\n");
  465. }
  466.  
  467. #undef LENGTH
  468. #undef NUMBER    
  469. #endif
  470.  
  471. static int pow5[21] [11] = {
  472.     {   1,   1,   0,   0,   0,   0,   0,   0,   0,   0,   0 },
  473.     {   1,   5,   0,   0,   0,   0,   0,   0,   0,   0,   0 },
  474.     {   1,  25,   0,   0,   0,   0,   0,   0,   0,   0,   0 },
  475.     {   1, 125,   0,   0,   0,   0,   0,   0,   0,   0,   0 },
  476.     {   2,   4, 113,   0,   0,   0,   0,   0,   0,   0,   0 },
  477.     {   2,  24,  53,   0,   0,   0,   0,   0,   0,   0,   0 },
  478.     {   2, 122,   9,   0,   0,   0,   0,   0,   0,   0,   0 },
  479.     {   3,   4,  98,  45,   0,   0,   0,   0,   0,   0,   0 },
  480.     {   3,  23, 107,  97,   0,   0,   0,   0,   0,   0,   0 },
  481.     {   3, 119,  26, 101,   0,   0,   0,   0,   0,   0,   0 },
  482.     {   4,   4,  84,   5, 121,   0,   0,   0,   0,   0,   0 },
  483.     {   4,  23,  36,  29,  93,   0,   0,   0,   0,   0,   0 },
  484.     {   4, 116,  53,  20,  81,   0,   0,   0,   0,   0,   0 },
  485.     {   5,   4,  70,   9, 103,  21,   0,   0,   0,   0,   0 },
  486.     {   5,  22,  94,  49,   3, 105,   0,   0,   0,   0,   0 },
  487.     {   5, 113,  87, 117,  19,  13,   0,   0,   0,   0,   0 },
  488.     {   6,   4,  56,  55,  73,  95,  65,   0,   0,   0,   0 },
  489.     {   6,  22,  26,  21, 112,  93,  69,   0,   0,   0,   0 },
  490.     {   6, 111,   2, 109,  51,  83,  89,   0,   0,   0,   0 },
  491.     {   7,   4,  43,  14,  35,   2,  34,  61,   0,   0,   0 },
  492.     {   7,  21,  87,  71,  47,  11,  44,  49,   0,   0,   0 }
  493. };
  494.  
  495. void pow_of5(int *v, int p)                                        /*;pow_of5*/
  496. /* This procedure is generated automatically 
  497.  * and therefore should not be modified.
  498.  */
  499. /* It has been modified for ANSI C */
  500. {
  501.     int i,*u;
  502.     if (p < 0 || p > 20) {
  503.         v[0] = v[1] = 1;
  504.         raise(SYSTEM_ERROR,"power of 5 too large");
  505.     }
  506.     u = pow5[p];
  507.     for (i=0;i<=u[0];i++) v[i] = u[i];
  508. }
  509.